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ABSTRACT 

Structural design problems can be considered to be optimization 
problems because a design is sought which is optimal by some criterion 
subject to limitations on size, behavior, or other aspects of the struc- 
ture. Under certain conditions, such problems may be solved by con- 
ventional mathematical programming techniques. The minimum weight 
design of a statically indeterminate three-bar truss is used to illustrate 
the application of the "sequential unconstrained minimization technique" 
of Fiacco and McCormick to the optimal design of trussed structures . A 
suggested computational procedure and computer program are included. 
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I. INTRODUCTION 



The desired conclusion of the structural design process is a structur- 
al configuration which performs its intended function efficiently. The 
method by which the final design is determined usually includes the 
establishment of an initial configuration from which the final design is 
evolved through a process of analysis and redesign. 

The electronic computer has permitted an acceleration of the design 
process by enabling the rapid analysis of the mathematical models assoc- 
iated with the structures under consideration. In order to more efficiently 
use this capability L. Schmit [1] has suggested a rational approach to the 
design process called "systematic structural synthesis." This approach 
involves the systematic evaluation and modification of a large number of 
trial designs until the optimal design is obtained. 

The basic problem is to design a structure which is optimal by some 
criterion subject to a set of requirements which specify acceptable limits 
on the behavior, size, or other aspects of the structure. Such problems 
are known to the operations analyst as mathematical programming problems, 
the general form of which is as follows . Find a vector X which 
minimizes F(X), 

subject to Hj (X) = 0, i = 1 ,2 , . . . ,m; 

Gj (X) £ 0, j = 1,2 ,p. 

T 

Here, X = (x 1( x 2 , . . . ,xj is an n-dimensional column vector; 

F(X) is the objective function; and the relations H t (X) = 0 and Gj(X) so 
represent the constraints on the problem . 
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If F(X) is convex in X and the constraint region is convex, i.e. the 
set of all solutions to Hj (X) = 0 , i = 1 , . . . ,m and Gj (X) = 0 , j = 1 , . . . , p 
forms a convex set, the problem is called a convex programming problem. 
If a problem can be classified as a convex programming problem then 
there are solution procedures which are guaranteed to find the global 
optimum. 

Several techniques which have been devised to solve this special 
case of the general programming problem are Rosen’s "gradient projection 
method," Zoutendijk's "method of feasible directions," and Kelly's 
"cutting-plane" method. These procedures are summarized by Hadley in 
Ref. 7. Another convex programming technique called the "method of 
alternate steps" was suggested by Schmit [l] as a means for converging 
on the optimal structural design. 

A recently developed convex programming technique which appears 
to be especially promising is the "sequential unconstrained minimization 
technique for convex programming with equality constraints" of Fiacco 
and McCormick [2] . By this method, the objective function and con- 
straints are dealt with simultaneously by combining them into a single 
function. This function is minimized for different values of an auxiliary 
parameter thereby generating a sequence of feasible solutions that con- 
verge to the optimal solution. 

If the convexity restrictions for the constraint region or the objective 
function are relaxed a convex programming solution technique may still be 
used but there is no guarantee that a solution to the problem is other than 
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a local minimum. For most practical problems, it is no less valuable to 



obtain information about local optima in the absence of a global solution. 

In general, structural design problems do not meet the convexity 
restrictions for convex programming problems . One such problem is the 
determination of the cross-sectional areas of the members of a statically 
indeterminant, coplanar, three -bar truss such that the total weight of the 
truss is minimized subject to constraints on stress and displacement. 
This structure is illustrated in Figure 1 . Procedures which can be used 
to solve this problem can also be applied to other trussed structures. 
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II. OBTECTIVE AND SCOPE 



The purpose of this study was to illustrate the application of the 
Fiacco-McCormick technique to the optimal design of trussed structures. 
Since minimum weight has usually been the major goal in the development 
of aircraft structural design techniques, it was used as the criterion for 
optimality . 

One example of the determination of the minimum weight design of a 
three-bar truss is included with results and the computer program. Also 
included is a brief discussion and outline of the solution procedure used 
to solve the problem. This procedure can only be applied to problems for 
which constraints can be formulated in completely analytical form. 
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HI. THE THREE-BAR TRUSS PROBLEM 



The three-bar truss problem of Schmit [l]is illustrated in the following 
sketch . y 




The function of the truss is to transmit the load Pj from the point 

of application at joint A to the support B. The deflections of the joint 

A in the x and y directions due to the load P, are u . and u . respectively. 

3 xj yj 

For a coplanar system of forces , one of the conditions for equili- 
brium of a rigid body in one plane is that the algebraic sum of the 
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components of all forces in each of two mutually perpendicular directions 
in the plane of the forces must separately vanish. Therefore, 

= p X j cos^ + p 2j cosjS 2 + p 3j cos/3 3 + Pj cosftj = 0 

and 



= p xj sinjS x + p 2j sin]3 2 + p 3j sinj3 3 - Pj sinotj = 0 
where p tJ is the force component in the direction of the ith structural 
member due to the jth load Pj . 

The stress in the ith member due to the jth force is defined as 
° ij = Pij /^i 

where A l is the cross-sectional area of the ith member. The equili- 
brium equations can then be written as 

A-lO-lj cos/^ + AjCr 2i cosjS 2 + A 3 ct 3j cos/S 3 + Pj cos&j = 0 (1) 

and 

A 1 o 1J sin Pi + A 2 a 2j sin|S 2 + A 3 o 3j sin£ 3 - Pj sinofj = 0 (2) 

Assuming small displacements, the change in length of the ith 
structural member due to the jth load is 



6 j, = - u . cos/Sj - u . sinpj 



xj 



yj 



(3) 



A relationship between stess, displacement, and temperature change can 
be established if is noted that may also be expressed as 

6,j = + o, AT, 1, . (4) 



where l t is the length of the ith member, oi t is the mean coefficient of 
thermal expansion of the ith member, E t is the modulus of elasticity of 
the ith member, and £T t is the temperature change imposed externally on 
the ith member. 
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If the temperature is assumed to be constant, equations (3) and (4) 

/ 

can be combined to yield 

Pij^i „ 

~ — — + u . cos p. + u . sinp, = 0 
Aj Ej xj ri uj 

and since = p^/Aj and 1 4 = N/sin/Sj , this equation becomes 
Ci, N . 

„ . — o~ + u . cos p. + u . sinp, = 0. (5) 

Ej sinPj xj 1 yj 1 

For purposes of this illustration, it is assumed that all geometric 

parameters of the truss are given except for the cross-sectional areas. 

It is further assumed that in addition to non-negativity restrictions for 

the Aj values, the behavioral limits for the structure are specified by 

upper and lower bounds on the stresses and displacements. Then, for 

the jth load condition, 

•^lj ^ , i = 1,2,3, 

L . £ u . £ U . , 

xj xj xj 

and 

L . £ u . £ U . . 
yj yj yj 

If pj denotes the specified density of the material in the ith struc- 
tural member then the weight, W, of the truss is 
W — Ay l x p x + A 2 I 2 Os A 3 1 3 p 3 
Since lj = I^/sin/Sj , this equation can be rewritten as 

■yy — A x N Py + A s N p 2 + A 3 N p 3 
sin£y sin£ 2 sin/5 3 

which is linear in A and is therefore convex. 
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The mathematical programming problem is then to fine A x , A 2 , A 3 , 
the cross-sectional areas of the structural members which minimize 

N* p 1 N* 02 N’Oa 

W = A, —77^+ As --- a- + As 3 - 



smf- 



sin^o 



1 sin/3 x 
subject to the constraints 

A]3i j cos + A^^ cos /? 2 + A 3 ct 3j cos £ 3 + cosc^ = 0 
A-lO^j sin/^ + A 2 ct 2j sin/S 2 + A 3 cr 3j sin/3 3 - Pj sinc^ = 0 
N 



^ E 1 sin/Sj. 



N 



2j E 2 sinjS, 



+ u .cosjSi + u .sinfi, = 0 
xj ^ xj 1 



+ u .cosjS 2 + u .sinjSo =0 
xj r * yj ^ 



N 



a3 J E 3 sin# 



+ u ,cosS 3 + u ,sin8 3 = 0 
, 3 xj P3 yj H3 



and 



Lj.^cr., * U., , L . * u . s U L , s u , s U ,,A,sO 
1 3 1 3 1 3 xj xj xj yj yj yj 

for i = 1,2,3 and j = 1 , 2 , . . . ,k where k denotes the number of distinct 
loading conditions. It is assumed that the loading conditions are not 
applied simultaneously. Note that inequalities such as L tJ s cr^ ^ Ujj 
can be rewritten as 

- Ljj ^ 0, 

-*11 + U u s 0. 

Now, in this problem statement, the variables to be determined are 
actually A x , A 2 ^ 3 , 0 ^ , i=l ,2,3, j = l , . . . ,k and u . , u . , j = 1 , . . . ,k . 

It is therefore apparent that the first two equality constraints are non- 
linear because of the cross products A } and hence the constraint 
region is non-convex. 
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IV. SOLUTION PROCEDURE 



Recall that the general form of the mathematical programming problem 
is to determine the vector X which 
minimizes F(X), 

subject to H t (X) = 0, i = 1,2, . . . ,m; 

Gj (X) *0, j = 1,2, . . . ,p, 

T T 

where X = (x x ,x 2 , . . . jXjj) 1 is an n-dimensional column vector. 

The basis of the Fiacco and McCormick technique is the P-function 
which is defined as follows for this general form. 

m p 

P(X,r) = F(X) + r L J (x) + A T (H^ (X) ) 2 . 

i=l 1 j=l 

The general solution procedure is then to seek minimum values of P(X,r) 

for r = r x ,r s , . . . ,r k where r x >r 2 . . . r k > 0 . Then 

lim X k = X 
k~> «° 

and 

lim P(X k ,r k ) = F 
k- 

where F is the minimum value of the objective function of the mathema- 
tical programming problem at X. 

The conditions for the determination of a global minimum for the P 
function are generally the same as for the convex programming problem. 
The objective function, F(X), must be a convex function. Recall that a 
linear function is convex although not strictly convex. Also the 
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constraint region must be convex (which is not true for this problem) 

and F , G x , . . . , G m , H x , . . . , H must be continuous . 

P 

Close inspection of the P function reveals that the term 
m 

r 1 Gj(X) 

i=l 

is a boundary repulsion term which keeps the minimum feasible solution 
of the P-function in the interior of the region defined by the inequality 
constraints. Naturally, the use of such a boundary repulsion term 
requires that solutions exist in the interior of the region defined by the 
inequality constraints . 

It is this feature of the Fiacco-McCormick technique which provides 
a significant improvement over other convex programming techniques. 

Other approaches, including Schmit's "method of alternate steps" require 
elaborate techniques to determine what action to take when the boundary 
is encountered during the minimization process. 

Another consequence of the boundary repulsion term which should be 
noted, however, is that the search for a feasible P-function minimum must 
be confined to the region defined by the inequality contraints . The 
reason for this restriction is apparent from the following example. Suppose 
that G x (X) = x x a 0. As x x -*0 from the negative side [G x (X)] -1 -* -<*>. 

The P-function can be minimized by direct analytical procedures for 
relatively trivial cases only. For other cases, some method of successive 
approximations must be used. One such method is that of gradient 
descent. 
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A. GRADIENT DESCENT 



Gradient methods are based on a Taylor's series approximation of 
the function to be minimized. 

1 . First-Order Gradient Method 

The first-order gradient method uses only the first-order terms of 
the Taylor's series expansion to obtain the equation, 

X 2 = X x - 0 7 P(Xj 

The procedure is then to move a distance proportional to 0 in the direc- 
tion of the gradient vector evaluated at the point X 1 . X 2 is the point at 
which the P-function is minimized along the gradient. 

2 . Second-Order Gradient Method 

A better approximation of the P-function is obtained by using the 
second-order terms of the Taylor's series expansion as well as the first 
to obtain the equation 



x 2 = x, - e 



t P(X x ) 



-1 



7P(X x ) 



dx t l Xj 

where the matrix || ♦ || is the Hessian matrix evaluated at X x and -1 
denotes the matrix inverse. The P-function is minimized along the vector 
given by the product of the inverse of the Hessian and the gradient of the 
function evaluated at X x . 

If the function to be minimized is not strictly convex, the Hessian 
matrix may be positive definite only in the vicinity of a local minimum. 
For this reason, the procedure used began with a search in the direction 
of positive 0 . If the P-function did not decrease in this direction before 
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encountering the boundary, the search was continued in the direction of 
negative & . This situation actually occurred in the example when the 
P-function was minimized at r=l . 

3 . Higher-Order Methods 

It would seem logical to consider the use of higher order terms to 
obtain an even better approximation of the function to be minimized. The 
effort required to program such a procedure, however, would be prohibitive 
for all but trivial problems. In fact, the labor required to calculate the 
second-order partial derivatives for the Hessian matrix may prohibit the 
use of the second-order gradient method for large problems . Experience 
with several small examples using both techniques indicates that the 
convergence properties of the second-order technique are far superior to 
those of the first-order technique. Fiacco and McCormick provide exam- 
ples of both techniques in Ref. 3 which clearly demonstrate the super- 
iority of the second-order method. 

B. ESTIMATION OF THE FINAL P-FUNCTION MINIMUM 

Fiacco and McCormick show in Ref. 2 that estimates of the solution 
of the P-function minimum for r=0 can be obtained by the use of poly- 
nomial approximations of x lr x 2 , . . . ,x^ and P as a function of r k ' 2 as 
follows. Expand each component of the X vector associated with the 
minimum of P(X,r) for a given value of r in a power series in r 1 ' 2 . To 
obtain a (p— 1 ) th -order polynomial approximation, drop the terms of the 
power series after the pth term. If the P-function has been minimized for 
r k , k=l , . . . , h where h ^ p , form the linear equations : 
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jzl 



_£lL 



Xl(r h-p } =a ° + --- +a ri( r h_n ) 



h-p 



1 + a i ( r u ) 

P-1 h-p 



x i( r J = a 0 +. . .+ aj_i(r ) 



j-1 

2 



_Rzi 



+ --- + a P -iK ) 



If it is assumed that r is reduced by a constant factor such that 

r k+ i = r k /C or r k = r x /C k-1 ( then the above equations constitute a system 

of p linear equations in p unknowns , a 0 , . . . , a . For the polynomial , 

P ^ 

2 P-1 

x t = a 0 + a x r + a^ + . . . + a^ ^r 

Xj at r=0 is given by a 0 . Hence, if these equations are solved for a 0 , 

a (p-1) th-order polynomial approximation is obtained for the value of x t 

corresponding to the solution to the mathematical programming problem. 

Estimates are obtained by this procedure for x, , . . . ,x as well as P at 

n 

r=0. Because of computer accuracy limitations, the order of the poly- 
nomial estimate will be restricted to three and will use the last four 



data points . 



C. CRITERIA FOR TERMINATION 

As r approaches zero, the value of the P-function minimum approaches 
the optimal solution. More rapid convergence to the optimal solutiom may 
be obtained by the extrapolation described above thereby permitting term- 
ination of the procedure with fewer iterations. The specific point at 
which the minimization process is to terminate depends, of course, upon 
the accuracy desired of the final solution estimate. 
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Fiacco and McCormick suggest several alternate criteria in Ref. 2 



based on the changes in the values of the variables or functions from one 
P-minimum to the next. In the case of the estimated value of the P- 
function minimum, the minimization process is terminated when 
| P * (r t— i) - P* (r t ) | < e 

where f is a small positive number and P* (r t ) is the polynomial estimate 
of the final solution after the P-function is minimized for the ith r value. 
For the three-bar truss problem, e was chosen to be 1x10" 5 . 

D. SUMMARY OF SOLUTION PROCEDURE 

The following procedures are those required to obtain an estimated 
value for the optimal solution to the mathematical programming problem . 

1. Form the P-function and obtain the analytical forms for the gradient 
vector and Hessian matrix. 

2 . Set p=4 for a third-order polynomial estimate of the values of P and 
x !,...,x for r=0 . 

3. Set r=l and h=l and select a starting point, X 0 , such that the in- 
equality constraints are satisfied. 

4. Determine the direction of steepest descent at the starting point. 
Search along this direction in the region defined by the inequality con- 
straints until the directional derivative of the P-function is less than the 
specified tolerance, i.e. near zero. 

5. Test the gradient magnitude. If less than a specified tolerance, go 
to step 8. 
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6. Call the present point the new starting point and go to step 4. 

7. If h s p calculate the (p— 1 ) th-order polynomial estimate of the com- 
ponents of the solution at r=0. 

8. If h ^ p+1 compare the last two estimates of the P function at r=0. 

If within the specified tolerance, terminate the routine. Otherwise, set 
r = r/C where C is the reduction factor and set h = h+1. Go to step 6. 
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V, THE EXAMPLE 



The example used to illustrate the application of the Fiacco- 
McCormick technique to the design of trussed structures is the problem 
described previously with two distinct load conditions. This yields a 
problem with 13 variables; 10 equality constraints, 4 of which are non- 
linear; and 23 inequality constraints including the three non-negativity 
constraints for the cross-sectional areas of the truss members . 

A. PROGRAMMING CONSIDERATIONS 

The computer program was written in Fortran IV language for the 
IBM 360 computer at the Naval Postgraduate School in Monterey, 
California. No attempt was made to optimize the program to reduce 
execution time. As it is written, the program required approximately 25 
seconds for execution. 

A simple bracketing procedure was used to locate the minimum of the 
P-function along the vector in the direction of steepest descent. The 
directional derivative was evaluated at the current starting point. A Q 
was then determined corresponding to a point at which the directional 
derivative was opposite in sign from that at the starting point and still 
within the region defined by the inequality constraints . The bracketing 
procedure was then commenced to reduce the magnitude of the directional 
derivative below the specified tolerance, DDTOL. 

Some experimentation was required to determine the necessary mag- 
nitude of DDTOL to ensure convergence of the gradient magnitude to within 
a tolerance of EPSI = lxl 0“ 3 . This was determined to be DDTOL = lxl 0~ 7 . 
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For the three-bar truss problem, a double-precision version of the 



program was required although several two-variable problems were 
solved without recourse to double precision. For small r values, the 
P-function becomes progressively more peaked in the vicinity of the 
minimum value. For the three-bar truss problem this peakedness is such 
that the change in ©required to obtain a bracket small enough to yield a 
directional derivative within tolerance is below the accuracy of the 
computer unless double precision is used. 

For programming purposes, the notation of the problem variables 

/ 1 / 2/ 3 / 11/ 21/ 31 / xl / yl / 1 2/ 22/ 3 2/ x2 ' y 2 was changed to 



F,x 1 ,x 3 , . . . ,x 13 respectively. 



B. INPUT DATA 

The input data for the problem was that of Schmit's Case 1 in Ref. 1. 

A cursory inspection will reveal that the values are highly unrealistic 

since the actual range of values of the modulus of elasticity is from 5 to 

30xl0 6 lbs ./in. 3 for metals. They are, however, sufficient for purposes 

of illustrating the solution technique. The input data is as follows. 

P 1 = 30 Pi = 1 jSj. = 135° E x = 1 

P 2 = 20 p a = 1 jS 2 = 90° E 2 = 1 

«i = 60° p 3 = 1 jS 3 = 45° E 3 = 1 

a 2 = 180° 

N = 1 

The constraints specified for the problem are 
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Al > ^2 1 ^3 ^ 0 / 



-15 s C7^ j j ^ 20, j 1,2, 

and 

-150 £ u .,u . ^ 200, j=l,2. 
x) yj 

T 

The starting point for the example was arbitrarily chosen to be X = 

T 

(2, 2, 2, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5) . The reduction factor for r was 1 0 at 
each iteration, i.e. r 1+1 =r t /10. 



C. RESULTS 

Values obtained for minimum P as a function of r as well as the third- 
order estimates of the P-function minimum for r=0 are shown in Table I. 
The value of the objective function, F, corresponding to the indicated r 
value; the total value of the remaining terms of the P function, denoted 
by P-F; and the values of the equality constraints, HjCX) are also shown. 

The value given under "Moves Required to Minimize" is the number 
of times a new direction of steepest descent was calculated prior to 
achieving a gradient magnitude of less than .001 for the current r value. 

The worth of the polynomial estimate is indicated by the fact that a 
final solution estimate of minimum P is obtained on the sixth iteration 
which differs by less than lxl 0~ 4 from the value obtained on the ninth 
and final iteration. 

By the ninth iteration, the value of the objective function, F, has 
converged to within 7x10“ 5 of the final solution estimate while the value 
of P(r 9 ) differs by 14x10" 5 from the final estimate. 
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Table I. Computer Solution for P-Function, Objective Function and Constraints 
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0.000000010 2 2.92253 2.92246 0.00007 -0.00000 -0.00000 0.00000 

3RD ORDER ESTIMATE 2.92239 
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Table II. Computer Solution for Cross-Sectional Areas, Stress, and Displacement. 
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Table II. (Continue) 
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One consequence of not obtaining an absolutely null gradient vector 
is revealed by the values of the equality constraints. Since convergence 
is obtained only to a gradient magnitude of .001, the equality constraints 
are not exactly satisfied for large r values. The values of the equality 
constraints are seen to decrease, however, at each iteration. In other 
words, the constraints get progressively "tighter" as r -* 0. The reason 
for this is the peakedness of the P-function in the vicinity of the minimum 
as r gets smaller and smaller. As P becomes more peaked, it is necessary 
to get nearer the minimum to satisfy the gradient magnitude tolerance 
for convergence . Thus the approximations become more accurate as 
r -» 0. A gradient magnitude tolerance of lxl 0 -4 failed to yield a 
change in the value of the constraints of sufficient magnitude to change 
the results although it did require a slightly higher execution time. 

Table II shows the values of the successive x 4 (r) , i=l , . . . , 13 as 
well as the final solution estimates. The final solution estimates changed 
by less than 2xl0“ 4 from the sixth to the ninth iteration. 

The values obtained by the method of Fiacco and McCormick are 
compared with those obtained by Schmit [1] using his "method of alter- 
nate steps" in Table III. 
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TABLE III 



Variable 


Fiacco & 
McCormick 


Alternate 

Steps 


h 


1.071 


1.072 


A 2 


.544 


.544 


A 3 


611 


.611 


^l 


19.877 


19.863 


*21 


20.000 


19.984 


°3 1 


.123 


1.206 


U xl 


19.755 


19.743 


u yi 


20.000 


-19.984 


CT 12 


-15.000 


-14.993 


*22 


5.000 


5.003 


a 32 


20.000 


19.996 


U x2 


-35.000 


-34.998 


U 9 

y2 


- 5.000 


- 5.003 


w 


2.922 


2.924 



The results obtained by the two methods are quite close although the 
minimum weight determined by the method of Fiacco and McCormick was 
slightly less. The primary reason for the difference was probably the 
termination criterion. Schmit solved the problem in 1960 on an IBM 653 
digital computer with severe storage capacity limitations so the accuracy 
of the solution in this study is unquestionably higher. Advances in 
hardware make any comparison of execution times meaningless. Schmit's 
results for c 31 is erroneous and is probably a typographical error in the 



reference . 
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VI. DISCUSSION 



A. ALTERNATE MINIMA 

A solution to the three-bar truss problem has been determined but 
the fact remains that it can be verified to be only a local minimum. Pro- 
cedures exist by which the "goodness" of the solution may be determined 
but the effort may be more than that required for solving the entire problem 
by a trial and error solution technique. One such procedure is simply to 
vary the starting point and search for alternate minima. 

The example problem was solved for several other starting points 
selected at random. These included 

X 0 = (1.0,1. 0,1. 0,5. 0,5.0,... ,5. 0) T 

= (0.5,0. 5,0. 5,0. 0,0. 0 ; ...,0.0) T 

T 

= (0.25,0.25,0.25,0.0,.. .,0.0) 

T 

= (0.25,0.25,0.25 ,-5.0, ... ,-5.0) 

= (0.25,0.25,0.25,-10.0, ... ,-10. 0) T 

The only change in value noted was a difference in the number of 
moves required to minimize the P-function for r=l . 

B. TWO -BAR TRUSS 

If the center bar of the truss is removed, the structure becomes a 
statically determinant two-bar truss in which the internal forces are 
uniquely determined. The minimum weight of such a structure subject to 
the constraints of the example is 3.04905. Thus, by adding a third 
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member, the total weight of the structure is reduced. It should be noted, 
however, that the behavior of the trusses will not be the same. 

Switsky [4] shows that the deflection of the two-bar truss will be 
less in this case than that of the three-bar truss. If the displacements 
of both trusses are required to be equal, given the same stress con- 
straints , the statically determinant two-bar truss will be the lighter of 
the two. 

C. N ON-ANALYTICAL CONSTRAINTS 

Structural problems often do not allow the nice problem statement 
form of the three-bar truss problem. The major difficulty arises from not 
being able to express some of the constraints in the completely analytical 
form required for the Fiacco-McCormick technique. Many problems 
require the use of empirically determined alignment charts or graphs to 
calculate changes in structural behavior as a result of the variation of 
design variables. For such problems, some variation of Schmit's 
"method of alternate steps" may be well suited. An example of such a 
problem is the determination of the minimum weight design of an integrally 
stiffened panel loaded axially under compression. Rosenbaum [5] des- 
cribes the problem and outlines a solution procedure based on the method 
of alternate steps. 
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VII. RECOMMENDATIONS FOR FURTHER STUDY 



Since only three design variables of the three-bar truss are allowed 
to vary, the example is actually a sub-optimization problem. This fact 
suggests that several possible extensions of the problem merit further 
consideration. 

Allowing more design variables such as the length of truss members 
to vary makes the problem even more nonconvex but the Fiacco-McCormick 
procedure may still be applied. The assumption of the availability of a 
continuous spectrum of materials would permit the determination of 
optimal material properties of the structure. 

The added restriction of a discrete choice of materials or material 
sizes available is a further realistic extension of the problem although a 
different solution procedure is required. 

Finally, the development of a systematic procedure for finding alter- 
nate minima is essential if one is to accept any result with confidence. 
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COMPUTER PROGRAM 



MAIN PROGRAM 

IMPLICIT FEAL^R (A-H,tl-Z) 

R E AL *8 Kll ,K12 ,K13,K21 ,K22 ,K23 ,K1A ,K15,K25 ,K31 ,KA2,K53 
$ ,K34,K3?,KVt,K4S,K54,K55.K24 
DIMENSION X( 13) ,XX(13, 12) ,XEST< 13,12) , X2 ( 13) ,SPR( 13, 13 
$ ) » COL ( 169) ,DELP( 13),X0LD( 13) ,LWK( 13) ,MWM 1 3) ,C ( 1 3 ) ,H ( 1 
$C , 12 ) 

D I ME NS I ON RR ( 12) ,P( 12 ) ,F( 12) » S I GMA (12),PEST(12),IA(12) 
DIMENSION XY ( 4 ) » RY ( 4 ) » PY ( 4 ) 

E QUI VALENCE! SPR <1 ),C0L(1 ) ) 

CCMMCM C1,C2,C3,K11,K12,K13,K21,K22,K23,K14,K24,K15,K2 
$5 * K3 1 » K42, K53 ,K34, K35, K44, KA5 » K54 » K55 
ODTOL - TOLERANCE FDR DIRECTIONAL DERIVATIVE. 

EPSI - TOLERANCE FOR GRADIENT MAGNITUDE 

DATA DDTOL/.CerOCOl/, EPS 1/ .00 1/ , DELTA/ 1. / 

STARTING X VALUES FOR SEARCH. 

DATA X/3*2. 0,10*5.0/ 

NV - NUMBER OF VARIABLES. 

NV= 1 3 
R R ( 1 ) = 1 . 

I R - NUMBER OF R VALUES FOR WHICH P FUNCTION IS TO BE 
MINIMIZED. 

I P=9 

J - NUMBER OF P VALUES FOR WHICH P HAS BEEN MINIMIZED. 

J = 0 

5 J = J + 1 

RR ( J ) - JTH P VALUE. 

R=RR ( J ) 

BEGIN SECOND ORDER GRADIENT SFARCH FOR MINIMUM OF P FOR 
CURRENT R VALUE. 

DDOLD - PREVIOUS DIRECTIONAL DERIVATIVE. 

DD0LD=C.C 

C GOLD - PREVIOUS GRADIENT MAGNITUDE. 

GCLD = 1 .C 

C XOLD - X VALUES CORRFSPONOI NG TO PREVIOUS GRADIENT 
C MAGNITUDE. 

DO 11 1=1, NV 
XOLD ( I )=r . c 
11 CONTINUE 

C II - COUNT OF MOVES IN SEARCH FOR MINIMUM OF P FOR CURRENT 
C R VALUE. 

11=0 

C EVALUATE HESSIAN MATRIX AT CURRENT STARTING POINT. 

C SPR - HESSIAN MATRIX. 

13 CALL MSOPA R ( X , R , S PR , NV ) 

C EVALUATE GRADIENT VECTOR AT CURRENT STARTING POINT. 

C DELP - GRADIENT OF P. 

CALL GRAD( X,R, DELP ) 

C INVERT HESSIAN MATRIX. 

CALL M INV( COL»NV»D»LWK»MWK) 

C STOP IF HESSIAN MATRIX IS SINGULAR. PICK NEW STARTING 
C POINT. 

IF(D.EO.F.C) WRITEI6, 17) 

17 FORMAT ( * *DET OF MATRIX OF 2ND ORDER PART I ALS = 0*') 
IF(O.FQ.r.C) STOP 

C PREMULTIPLY INVERSE OF HESSIAN BY GRADIENT VECTOR. 

CALL MTXMULI SPR, DELP, NV,C) 

C GMAG - MAGNITUDE OF GRADIENT. 

GMAG=C .C 
DO 23 1=1, MV 

GMAG=GMAG+(DELP( I )*0ELP( I ) ) 

23 CONTINUE 

GMAG=D SQRT(GMAG) 

C TEST MAGNITUDE OF GRADIENT. IF BELOW TOLERANCE, P IS 
C MINIMIZED FOR CURRENT R. 

IFIGMAG.GT .EPSI ) GO TO 28 
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C PSTAR - P MINIMUM FOR CURRENT R VALUE. 

PSTAR=PVALUE (X ,R) 

GO TO 85 

C IF GR A C I ENT MAGNITUDE IS SMALL (LESS THAN DELTA) AND 
C GMAG INCREASES, CONVERGENCE IS UNCERTAIN. TERMINATE 
C SEARCH WITH SMALLEST GRADIENT MAGNITUDE OBTAINED. 

28 IF (GMAG. LT .DELTA. AND. GMAG. GE. GOLD ) GO TO 30 
GO TO 38 

30 WRITE(6,3I ) GOLD 

31 FORMAT ( • ^MAGNITUDE OF GRADIENT WILL NOT CONVERGE TO • 
i'WITHIN TOLERANCE. LAST GMAG OBTAINED WAS ' F10 . 5» * ' / ) 

DO 35 1=1, NV 
X( I ) =XOLO ( I ) 

35 CGNTINUE 

PSTAR=PVALUF (X ,R) 

GO TO 85 

C SAVE PRFSENT GRADIENT MAGNITUDE AND LOCATION FOR NEXT 
C CONVERGENCE TEST. 

38 GOLD=GMAG 
DO 41 1 = 1, NV 
XOLD ( I ) = X ( I) 

41 CONTINUE 

C EVALUATE DIRECTIONAL DERIVATIVE AT STARTING POINT. 

C OD - DIRECTIONAL DERIVATIVE. 

DD=D I R DEE ( NV , DFLP , C ) 

C IF DIP ECT I ON AL DERIVATIVE IS WITHIN TOLERANCE BUT GRADIENT 
C IS NOT, CONVERGENCE IS IMPOSSIBLE WITH PRESENT TOLERANCES. 
IF(DABS( DDJ.GT .DDTOL) GO TO 48 
WRITE(6,45 ) DD , R , I I 

45 FORMAT ( • * EPSI TOO SMALL OR DDTOL TOO LARGE. SEARCH W' 
t'ILL NCT CONVERGE. LAST DD OBTAINED = • F 1 2. 8 • . R = • F 10 . 3 
* • . I = ' I4**» ) 

STOP 

C 

C BEGIN GRADIENT SEARCH FOR MINIMUM IN CURRENT DIRECTION OF 
C STEEPFST DESCENT. 

C 

C UNITIZE DIRECTIONAL DERIVATIVE. KEEP SIGN. EVALUATE P 
C P UNCT I ON AT START POINT. 

48 TEST=DD/DABS (DD) 

PV1 = PVALUF (X ,R ) 

C INITIAL THETA VALUE IS POSITIVE. 

T1=0.C 
STFP=. 5 
52 T=T1+STEP 

C STEP IN DIRECTION OF STEEPEST DESCENT TO GET X2. 

CALL XVALUE( X2 ,X,NV,C,T) 

C TEST X2 FOR FEASIBILITY WITH RESPECT TO INEQUALITY 
C CONSTRAINTS. 

CALL CHECK ( X2, NV, RESULT) 

C IF RESULT IS ONE, X2 IS FEASIBLE. IF X2 IS NOT FEASIBLE, 

C REDUCE THETA VALUE. 

IF ( RESULT . EQ . 1 .' ) GO TO 58 
STEP=STEP/2. 

GO TO 52 
58 T 1=T 

C EVALUATE P AFTER STEP. 

PV 2=PVALUE ( X2, R> 

C EVALUATE 0 IR EC T I ON AL DERIVATIVE AT X2. 

CALL GRAD( X2,R,DELP) 

DD=DIRDEP ( NV , DELP ,C ) 

C TEST FOR MINIMUM AT X2. 

IF(DABS(DD).LE. DDTOL) GO TO 82 
C TEST SIGN OF DIRECTIONAL DERIVATIVE AT X2. IF OPPOSITE 
C FRCM THAT AT STARTING POINT, BRACKET HAS BEEN ESTABLISHED. 
TR Y=DD /DAB S( DD ) 

IF (TRY+TEST) 65,68,65 

C IF NCT OPPOSITE IN SIGN BUT P INCREASES, SEARCH FOR 
C MINIMUM IN DIRECTION OF NEGATIVE THETA FROM START POINT. 

65 IF(PV2.LT.PV1) GO TO 52 
STEP=-2. 

GO TO 52 
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C INITIAL BRACKET IS ESTABLISHED. REDUCE STEP SIZE BY 1/2. 

C REVERSE DIRECTION OF SEARCH WITH REDUCFD STEP SIZE. 

68 TE ST=T RY 
STEP=T 

70 STEP=-STEP/2. 

71 T=T+STEP 

CALL XVALUEIX2,X,NV,C,T) 

CALL GRADI X2 ,R ,DELP) 

DD=DIRDER< NV,DELP,C) 

C TEST FOP MINIMUM AFTER EACH STEP. 

IF(DABSIDD).LE.DDTOL) GO TO 82 
C TERMINATE SEARCH IF DIRECTIONAL DERIVATIVE REPEATS. 

IF (DO. EQ.DDOLD ) GO TO 82 
C SAVE CURRENT VALUE OF DD. 

DDCLD= DD 

C TEST FOR ESTABLISHMENT OF REDUCED BRACKET. 

TP Y=DD /DAB S( DD) 

IF I TRY +TEST ) 71,80,71 
80 TE ST = TRY 
GO TO 7C 

C CALCULATE VALUE OF NEW STARTING POINT. BEGIN SEARCH IN 
C NEW DIRECTION OF STEEPEST DESCENT. 

82 CALL XVAl UE(X, X,NV,C,T) 

11=11+1 
GO TO 13 
85 I A ( J )= 1 1 
P( J) =P STAR 

C EVALUATE ORIGINAL OBJECTIVE FUNCTION. 
F(J)=XI1)*C1+XI2)*C2+XI3)*C3 
SI GMA ( J)=P( J)-F( J) 

C EVALUATE EQUALITY CONSTRAINTS. 

H( 1, J)=X{ 1 )*X( 4)*Kll+X( 2)*X( 5) *K12+X( 3)*Xl 6) *K13-K14 
H{ 2, J) =X( 1 )*X< 4)*K21+X{ 2>*XI 5)*K22+X<3 ) *X I 6 ) *K23-K24 
H( ?, J)=K31*X< 4)+K34*X< 7)+K3 5*X< 8) 

H( A, J)=K42*X (5 )+K44*Xl7)+K45*X( 8) 

H ( 5, J)=K53*X< 6)+K54*X( 7)+K55*X< 8) 

HI 6, J >=X{ 1 )*X< 9 )*K 1 1+X { 2 )*X( 10 )*K1 2+XI 3) *X I 1 1 ) *K13-K15 
HI 7, J) =XI 1 )*X< 9)*K21+XI2)*X< 10)*K22+X( 3)*X 11 1 )*K23-K25 
HI 8, J)=K31*X<9 )+K34*X< 12 ) +K3 5* X I 1 3 ) 

H I 9 , J ) = K42 *X I 1 0 ) +K44*X I 12)+K45*X( 13) 

HI 10,J)=K53*X( 11)+K54*XI 12)+K55*X< 13) 

DO 92 K= 1 , NV 
XX IK, J )=Xl K) 

92 CONTINUE 

IFIJ.LT.4) GO TO 108 

BEGIN THIRD ORDER POLYNOMIAL CURVE FIT TO ESTIMATE P AND X 
VALUES FOP R=C . PY,RY,AND XX CONTAIN LAST FOUR VALUES OF 
P , R , AND X IN F CRM REQUIRED FOR SUBROUTINE POLY. 

u=c 

J3=J-3 

DO 10O K= J 3, J 
11=11+1 
PY I I 1 ) =P I K ) 

RYI I 1 ) =RR I K) 

100 CONTINUE 

CALL PCLYI PY,RY,4, AO) 

C PEST - ESTIMATE OF P AT R=0. 

PEST I J ) = AO 
DO 105 1=1, NV 
I 1=0 
J3=J-3 

DO 1C2 K= J 3, J 
I 1 = 1 1+1 

XYI II )=XX( I , K ) 

102 CONTINUE 

CALL POLYI XY,RY,4, AO) 

C XE ST - ESTIMATE OF X AT R = 0. 

XESTI I , J) = AO 
105 CONTINUE 

IFIJ.EQ.IR) GO TO 122 
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108 J1=J+1 

C REDUCE R BY ARBITRARY CONSTANT. 

RR ( J 1 ) =RR ( J ) / 1 0 . 

GO TO 5 

122 WR ITE ( 6, 12 3 ) 

123 FORMAT (• 1 • 44X • MOVE S * /A 3X • REQU I R ED ' /A6 X • TO • /23X'LINE*8X 
$• R *7X* MINI MI ZE • AX* P* 8X* F* 7X* P-F*8X*H1 ' 8X'H2' 8X*H3U 

DO 133 J=1,IP 

WPITE(6,128)J,RR(J), I A ( J ) , P C J ) , F ( J ) , S I GM A ( J I , ( H( I* Jit I 
$=1,3) 

128 FORMAT (*C*22X,I3,3X,F12.9,3X,I3,3X,F8.5,F9.5,F1?.5, 
$F1C.5»F1C.5»F10.5) 

IFIJ.LT. A) GO TO 133 
WRITE(6,132) PEST(J) 

132 FORMAT ( * *29X'3RD ORDER ESTIMATE *F8.5) 

133 CONTINUE 
WRITE16, 135) 

135 FORMAT (• 1 • 23 X • L I NF • 7X • R • 1 A X ' HA • 8X • H5 • 8 X • H6 • 3X ' H7 • 8X • H8 
$8X *H9 ' 8X *H ID ' ) 

WRITE (6,138) (J,RR(J)» ( H ( I , J) , I =A , 1 " ) , J = 1 , I R ) 

1 3 8 FORMAT! • C • 23X , 1 3 , F 1 A. 9 , 7F 1 0 . 6 > 

WR ITE ( 6* 1 AC ) 

1 AO FORMAT (• 1' 23X ' L I NE • 7X • R • 1 1 X * X 1 • 8X ' X2 * 8X ' X 3 * 3X • XA • 8X • X 5 
$8X *X6‘ 8X »X7' ) 

DO 1 AR J=1 , IR 

WRITE (6,1 A A) J,RR( J) , (XX! I , J) , 1=1,7) 

IAA FORMAT! 'C • 23X, I 3, F 1A. 9, 7F 10. 5) 

IF(J.LT.A) GO TO 1A8 

WRITE! 6, 1A7) ( XE S T ( I , J ) , I = 1 ,7 ) 

1 A 7 FORMAT!* »,28X*3RD ORDER E S T • F 9. 5 , 7F 1 0 . 5 ) 

1A8 CONTINUE 

WRITE! 6, 150 

150 F0RMAT(»1*28X , LINE'7X , R'11X , X8'3X , X9‘8X , X10*7X , X11*7X 
$ ' X 12 ' 7X • XI 3' ) 

DO 158 J = 1 , IR 

WRITE (6, 15A) J,RR( J ) , ( XX! I , J) , 1=8, 13) 

1 5 A FORMAT ( , 0 , 23X,I3»F1A.9,6F10.5) 

IF(J.LT.A) GO TO 158 

WRITE! 6, 15 7) (XEST! I, J) , 1 = 8, 13) 

157 FORMAT!* *,33X*3RD ORDER E ST • F9 . 5 , 6F 1 r . 5 ) 

158 CONTINUE 
STOP 
END 



SUBPROGRAMS 



SUBROUTINE POL Y ( X , RR , K , AO ) 

THIS SUBROUTINE SOLVES THE LINEAR EQUATIONS FOR THE 
POLYNOMIAL ESTIMATE OF THE VARIABLE CONTAINED IN 
ARGUMENT X. AO IS THF ESTIMATED VALUE OF THE VARIABLE 
AT R = r . 

IMPLICIT P EAL*8( A-H,0-Z) 

DIMENSION X(A),RR(A).A(A,A),L(16),M!16),C(16) 
EQUIVALENCE! A! 1) ,C(1) ) 

DO 1 J = 1 , K 
DO 1 1=1, K 

A! I, J) =(DSQRT(RR( I )) )**( J-l) 

1 CONTINUE 

CALL M INV ( C, K, D, L , M) 

AO =( .0 

OC 2 J 1= 1 , K 

AO = AO+ ( A ( 1 ,J1 ) )*X( J1 ) 

2 CONTINUE 
RFTURN 
END 



DOUBLE PRECISION FUNCTION D I RDER ( N V , DELP , C ) 

C THIS SUBPROGRAM CALCULATES THE VAL"E OF THE DIRECTIONAL 
C DERIVATIVE. 
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IMPLICIT REAL*8(A-H»0~Z) 

RE AL*8 DFLP(NV) ,C ( NV ) ,MAGSQ 

MAGSQ=O.C 

DO 1 1=1, NV 

MAGSQ=MAGSQ+(C(I)*C(I)) 

1 CONTINUE 

IF (MAGSQ. EQ.O .0 ) GO TO 3 
UV=1/DSQRT (MAGSQ) 

D IRDER=C .0 
DO 2 1=1, NV 

DIRDER=DIRDER+ <DELP< I )*C( I )*UV) 

2 CONTINUE 
RETURN 

3 DIRDER =C.O 
RETURN 

END 



SUBROUTINE MTX MUL ( SPR , DEL P , NV , C ) 

C THIS SUBROUTINE PREMULTIPLIES THE VFCTOR UELP BY THE 
C MATRIX SPR. 

IMPLICIT P EAL*8( A-H,0-Z) 

DIMENSION SPR ( NV » NV ) ,DELP(NV),C(NV) 

DO 1 I = 1 , N V 
C( I)=O.C 
DO 1 J = 1 »N V 

C( I)=C( I )+(SPR( I , J)*DELP< J) ) 

1 CONTINUE 
RETURN 
END 



SUBROUTINE XVALUE ( X2 . X 1 , J , C , T ) 

C THIS SUBROUTINE CALCULATES THE POINT X2 A DISTANCE 
C PROPORTIONAL TO T FROM Xl IN THE DIRECTION OF 
C STEEPEST DESCENT. 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION X2( J),X1 (J) ,C(J, 1) 

DO 1 I = 1 ,J 
X 2 ( I ) = XI (I )- <T*C( I, 1 ) ) 

1 CONTINUE 
RETURN 
END 



SUBROUTINE CHE CK ( X , NOV AR , RESULT ) 

C THIS SUBROUTINE CONTAINS THE INEQUALITY CONSTRAINTS. 

C IF THE POINT X SATISFIES THESE CONSTRAINTS, RESULT=+1. 
IMPLICIT REAL*8( A-H,0-Z) 

REAL *8 X< 13),L4,L5,L6,L7,L8,L9,L10,L11 ,L12 ,L13,K11 ,K12 
$ , K13,K21 ,K22,K23,K14,K15,K25,K31,K42,K53,K34,K35,K44, 
$K45,K54,K55,K24 

COMMON Cl, C2,C3,K1 1,K12,K 13 , K 2 1 , K2 2 , K2 3 , K 1 4 , K24,K15, 
$K25,K31,K42,K53,K34,K35,K44, K45,K54,K55,U4,U5,U6,U7,U8 
$,L4,L5,L6, L7, L 8 

EQUIVALENCE (L4,L9I,(L5,L1G),(L6,L11),(L7,L12),(L8,L13 
$) , (U4.U9) , (IJ5.U10) , (U6 ,U11 ) , (U7,U12» , (UR,U 13) 

I F ( X ( 1 ) .Lt . 3 .0 ) GO TO 1 
IF(X(2).LE.9.0) GO TO 1 
IF(X(3).LE .0.0) GO TO 1 
IF(X(4).LE.L4.0R.X(4) .GE.U4) GO TO 1 
IF(X(5).LE.L5.0R.X(5).GE.U5) GO TO 1 
IF(X(6).LE .L6.0R. X(6) .GE.U6) GO TO l 
IF(X(7).LE.L7.0R.X(7).GE.U7) GO TO 1 
IF(X(8).LE.L8.0R.X(8).GE.U8) GO TO 1 
IF (X( < 3).LE.L9.0R.X(9) .GE.U9) GO TO 1 
IF(X(1C) .LE.LIO.OR.XI 1O.GE.U10) GO TO 1 
IF(X(11).LE.L11.0R.X(11).GE.U11) GOTO 1 
I F < X < 12) .LE.L12.0R.XI 12).GE.U12) GO TO 1 
IF(X (13) .LE.L13.0R.X( 13) .GE.U13) GO TO 1 
RE SULT = 1.0 
RETURN 
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I RE SULT =—1.0 
RETURN 
END 

DOUBLE PRECISION FUNCTION PVALUE(X,R) 

C THIS SUBPROGRAM CONTAINS THE ANALYTICAL FORM OF THE 
C P-FUNCTION • THE VALUE OF P AT X l S CALCULATED. 

IMPLICIT REAL*8< A-H.O-Z) 

REAL*8 X(13) |L4,L5»L6, L7,L8,L9,L10,L11 »L12 ,L13*K11,K12 
$, K13,K21 ,K22tK23,K14,K15,K25,K31,K42,K53 ,K34,K35,K44, 
$K*5»K54,K55,K24 

COMMON C1,C2,C3,K11,K12,K13,K21,K22,K23,K14,K24,K15, 
$K25, K3I,K42,K53,K34,K?5,K44, K45 , K 54, K 5 5, U4 , U5 , U6 , U7 , U8 
* t L4.L5,Lfc, L7.L8 

EQUIVALENCE (L4 f L9),(L5,L10),(L6,Lll),(L7,L12),<L8,L13 
$) , (U4,U9) , <U5,U10) , (U6,U11) , <U7,U12) , (U8,U13) 

Hl = Xm*X( 4) *K11+X(2> *X(5) *K12+X<3 )*X (6 ) *K 13-K14 
H2=X( 1 )*X< 4) *K21 + X( 2) *X ( 5 ) *K 22+ X ( 3 ) * X( 6 ) *K2 3-K24 
H3=K3 1 *X (4 ) +K34*X ( 7)+K35*X( 8) 

H4 = K42*X ( 5 ) + K44*X( 7)+K45*X( 8) 

H5=K53*X(6 )+K54*X( 7)+K55*X( 8) 

H6 = X(1I*X( 9) *K11+X(2)*X( 10 ) *Kl 2 +X ( 3 ) *X ( 11 ) *K13-K15 
H7=X ( 1 )*X( 91 *K21+X(2)*X( 10 ) *K2 2+X ( 3 ) * X ( 1 II *K23-K2 5 
H8=K31*X (9 ) +K34*X ( 12) +K35*X( 13) 

H9=K42*X( 1 0) ♦K44*X { 12 ) +K45*X(13) 

H1C=K53*X( 11 ) +K54* X ( 1 2 ) +K55*X( 1 3) 

G 1 =1 . / X ( 1 ) 

G2=l./X( 2) 

G3=l./X( 3) 

G4=l./ ( X ( 4 )-L4 > 

G5=l • / ( XI 4 )-U4 ) 

Gfc-1./ ( X < 5 )-L 5 ) 

G7=l . / ( X ( 5 ) -U5 ) 

G8=1./(X(6)-L6) 

G9 = l ./ (X I 6 )-U6 ) 

G10-1. /(X( 7) -L7) 

G11=1./(X( 7I-U7) 

G1 2*1 • / ( X ( 8I-L8) 

G 1 3= 1 . / ( X ( 8) -U 8 1 
G14=1./(X( 9J-L9) 

G15=1./(X( 9I-U9) 

G 1 6= 1 . / ( X( 1 J1-L10) 

G17*1./(X( 10 )-U10 ) 

G 1 8= 1 . /<X( 11 l-Lll ) 

G19*1./(X( lll-Ull) 

G20=1./(X< 12 I —LIP » 

G2 1= 1. /<X< 12 ) -U12 ) 

G2 2= 1 . / ( X { 13 ) -L 13 I 
G2 3= 1 . / ( X ( 13 ) -U1 3 ) 

Y 1= 1 . /DSQP T ( R ) 

EQ=H 1**2 +H2* *2 +H3 **2+H4**2+H 5**2+H6** 2+H7* *2 + H8** 2 + 
SH9**2+H10**2 

UN=G 1+G2+G 3+G4-G5+G6-G7+G8-G9+G 10-G 11+G12-G1 3 + G14-G1 5 ♦ 
$G16-G17+G1P-G19+G20-G21+G22-G23 
PV ALUE*C 1* X ( 1 ) +C2*X(2)+C3*X(3)+Y 1*EQ+R*UN 
RETURN 
END 



SUBROUTINE GRAD( X* R.D) 

C THIS SUBPROGRAM CONTAINS THE ANALYTICAL FORM OF THE 
C GRADIENT VECTOR. THE VECTOR IS EVALUATED AT X. 

IMPLICIT REAL*8( A-H.O-Z) 

REAL*8 D( 1 3 ) 

PEAL*8 X(13),L4,L5 f L6,L7,L8,L9,LlC,Lll,L12,L13tKll,K12 
$, K13,K21,K22,K23,K14,K15,K25,K31,K42,K53,K34,K35 , K44 , 
$K45,K54,K55,K24 

COMMON CltC2»C3,Kll,K12tK13fK21fK22fK23.K14,K24,K15, 
$K25,K31,K42,K53,K34,K35,K44,K45,K54,K55,U4 ,U5, U6,U7,U8 
$,L4,L5,L6, L7.L8 

EQUIVALENCE (L4,L9).(L5,L1C).<L6,L11>,(L7,L12),(L8,L13 
$) , (U4,U9) , (U5» U10 ) , (U6,U11 ) , (U7.U12) , (U8,U13) 

DIMENSION R0(3) ,E(3) , B A < 3 ) , AA < 2 > , P { 2 ) 
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H1=X( 1 )*X( 4)*K 11+ XI 2) *X( 5) *K12+X< 3 ) * X < 6 ) * K 1 3-K14 
H2=X ( 1 )*X< 4)*K21+X(2)*X(5)*K22+X< 3)*X( 6)*K23-K24 
H3=K31*X(4)+K34*X(7)+K35*X{8) 

H4=K42*X( 5 )+K44*X( 7)+K45*X( 8) 

H5 =K53*X ( 6 ) +K5 4*X ( 7)+K55*X( 8) 

H6=X( 1 ) *X( 9)*K11+X(2)*X( 10 ) *K1 2 + X ( 3 ) *X ( 1 1 ) *K 13-K1 5 
H7=X(1)*X<9) *K 21+ X ( 2 ) *X( 1 C ) *K2 2+X ( 3 ) * X ( 1 1 ) *K2 3-K2 5 
H8=K31*X(9)+K34*X(12)+K35*X< 13) 

H9=K42*X{ 1C)+K44*X< 12 )+K45*X( 13) 

Hi 0=K5 3*X ( 11 )+K54*X< 12)+K55*X( 13) 

Gl=l./X( 1 ) 

G2=l./X( 2) 

G3= 1 . / X ( 3 ) 

G4=1./(X(4)-L4) 

G5=l • / ( X { 4 )-U4 ) 

G6=1./(X(5 )-L5 ) 

G7=l ./ ( X ( 5 )-U5 ) 

G8=1./(X(6)-L6) 

G9=l ./ (X (6 )-U6 ) 

G10=1./(X< 7) -L 7) 

G 1 1= 1 . / ( X ( 7)-U7) 

G 1 2=1 • / ( X ( 8J-L8) 

G 1 3= 1 . / ( X( 8) -U 8 ) 

G14=1./(X(9)-L9) 

G 1 5= 1 . / ( X < 9 ) -U9 ) 

Glfc= 1. /( X( 13 ) - L 1 0 ) 

G1 7= 1 • / ( X ( 10 )-U10 ) 

G 1 8= 1 . /(X( 1D-L11) 

G19=1./(X( 1D-U11) 

G2 C= 1 • / { X ( 12 )-L12 ) 

G21=l. /( X( 12 ) - U12 ) 

G2 2= 1 . / ( X ( 13 )-Ll 3 ) 

G23=l# /C X ( 13 )-U13 ) 

Y2=2 • /DSQP T( R ) 

D( 1)=C1+Y2*(X (4)*( K11*H1+K21*H2)+X(9)*(K11*H6+K21*H7) ) 
$-R*(Gl**2) 

D ( 2 ) =C 2+Y 2 *( X( 5)*< K12*Hl+K22*H2) + X( 10)*( K12*H6+K22*H7) 
^ j _ ( G2 a * c3 ^2 ) 

D( 3)=C 3+Y2*< X ( 6)*( K13*H1+K23*H2 ) + X ( 1 1 ) *< K1 3*H6+K23*H7 ) 
$ ) -P* { G 3**2 ) 

D< 4)=Y2*(X (1 )*(K11*H1+K21*H2)+K31*H3 > -R* ( G4**2-G5**2 ) 

D ( 5)=Y2*(X (2)* (K12*H1+K22*H2)+K42*H4) -R* < G6**2-G7**2 ) 
D( 6)=Y2*(X (3 )*<K13*H1+K23*H2 )+K53*H5)-R*(G8**2-G9**2) 
0( 7)=Y2*(K34*H3+K44*H4+K54*H5 )-R*(G 10 **2-G 11**2 ) 

D( 8)=Y2*(K35*H3+K45*H4+K55*H5)-R*(G12**2-G13**2) 

0 ( 9 ) =Y 2* ( X (1 )*<K1 1*H6+K21*H7)+K31*HS )-R*( G 14**2-G 15**2 
$) 

D ( 1C )=Y2*( X( 2) *(K12*H6+K22*H7)+K42*H9)-R*( G16**2- 
$G17**2 ) 

D ( 1 1 ) = Y2*( X(3)*(K13*H6+K23*H7)+K53*H10 )-R*(G18**2- 
$G19**2 ) 

0(12 ) = Y2* ( K34*H8+K44*H9+K54*H10 ) -R* ( G20**2-G2 1**2 ) 

D ( 13 ) = Y2*( K35*H8+K45*H9+K55*H10)-R*(G22**2-G23**2 ) 

RETURN 

END 



SUBROUTINE MSOPAR < X.R , S, NV ) 

C THIS SUBROUTINE CONTAINS THE ANALYTICAL FORM OF THE 
C HESSIAN MATRIX. THE MATRIX IS EVALUATED AT X. 

C PROBLEM CONSTANTS ARE READ INTO COMMON IN THIS ROUTINE. 

IMPLICIT REAL* 8( A-HtO-Z) 

RE AL*8 S ( 1 3» 13) 

RE AL*8 N 

REAL *8 X< 13) ,L4,L5.L6,L7,L8,L9,L10,L11 1 L12 ,L13,K11 ,K12 
$»K13 tK21,K22,K23,Kl4,K15,K25,K31,K42,K53,K34,K35,K44, 
$K45,K54,K55,K24 

COMMON C1»C2,C3,K11,K12»K13»K21 , K2 2 , K23 , K1 4, K24, K 1 5 , 
$K25 f K31,K42,K5 3 ♦ K 34 » K35 , K44 , K45 , K54, K 55 , U4 t U5, U6 , U7, U8 
$,L4 t L5,L6,L7,L8 

EQUIVALENCE ( L 4. L 9 ) , < L 5 , L 1 0 ) . < L6, L 1 1 ) , ( L 7 , L12 ) , ( L8 , L 1 3 
t) ,(U4,U9), (U5 »U10 ) »(U6»U11)» (U7»U12)» ( U8 » U 13 ) 
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01 HENS ION R0I3) ,E13) ,BAI3) ,AAI2),P(2) 

OAT A RO ( 1 ) ,Rl)(2)»R0ni,EI l ) » E l 2 ) » 1 1 3 ) , N/7* 1. 0/ 

DATA L/C/ 

L = L+ 1 

IF(L.GT.l) GO TO 2 

114=20.0 

U5=20 . 0 

U6=20.C 

U7=2CC. 

U8=20C . 

L4=- 1 5 . 

L 5=- 15. 

LG=- 15. 

L7=-2CC . 

L8=- 200. 

PI=3. 1415926536 
BAIl ) = .75*PI 
BAI2) = .5C*PI 
BA 13 )= .25*PI 
AAIl)=PI/3. 

AA I 2 ) = PI 
PI 1 ) = 3C . 

P ( 2) =2C. 

C 1=N*R0I 1) /OS INI BA ( 1 ) ) 

C2 =N*P 01 2I/0SINIBAI2) ) 

C3=N*R013> /OS I N( BA ( 3 ) ) 

K 1 1=0C0S I B AI 1) ) 

K12 = 0C0S(BAI2) ) 

K 1 3=DC0S I B AI 3) ) 

K21=DS INI B AI 1 ) ) 

K22 = 0SINIRAI2) ) 

K23=DS IN I 3 AI 3 ) ) 

K14=-P(1)*DC0S(AA( 1 ) ) 

K24=PI 1)*DSINI AAI 1 > ) 

K15=-PI2)*DC0S(AAI 2) ) 

K25=PI2)*DSIN( AA I 2 ) ) 

K31=N/(E(1 )*DSIN(BA(i) ) ) 

K4 2=N/ I E I 2 )#DS INI BAI2) ) ) 

K5 3=N/ IE 13 ) *DS I N I BA I 3 ) ) ) 

K34=K 1 1 
K35=K2 1 
K44=K12 
K45=K2 2 
K 5 4= K 1 3 
K5 5= K2 3 

2 Y2=2 . /OS ORT I R ) 

H 1 = X I 1 )*X< 4)*K11+X12)*X( 5) *K12+X(3 ) *X < 6 ) *K 13-K14 
H2=X( 1 )*X I 4)*K 2 1 + X I 2) *Xt 5) *K22 + XI 3 >*Xl 6)*K23-K24 
H3=K31*X(4)+K34*X(7)+K35*X(3) 

H4=K42*XI 5)+K44*Xl 7)+K45*X( 3) 

H5=K53*X(6 )+K54*X( 7)+K55*X( 3) 

H6 = X 1 1 )*X(9)*K11 + X(2)*XI10)*K12+X(3)*X(11)*K13-K15 
H7=X I 1 ) *X( 9) *K21 + X< 2) *X( 1G ) *K2 2 + X I 3 ) * X I 1 1 ) *K23-K25 
H8=K31*X 19 )+K34*X I 12)+K35*X{ 13) 

H9=K42*X( 1 0 ) +K44* X(12)+K45*X(13) 

HI G= K5 3*X ( 11 )«-K54*X( 12 ) +K55*X( 13) 

Gl=l./X{ 1 ) 

G2=i./X(2) 

G3 = 1 . / X I 3 ) 

G4 = l . / I X I 4 ) -L4 ) 

G5 = l . / I XI 4 )-U4 ) 

G6 = l . / I X I 5 )-L 5 ) 

G7=l •/( XI 5 )-U5 ) 

G8=1./(X(6)-L6) 

G9 = l ./ I X 1 6 )-U6 ) 

G 1 0= 1 . / 1 XI 7 ) -L 7 ) 

G 1 1= 1 • / 1 X I 7J-U7) 

G12=1./(X( 8 ) -L 8 ) 

G 1 3= 1. /( XI 8 ) -U 8 ) 

G14=1./(XI 9)-L9) 

G15=1./(X(9)-U9) 

G16=l./( XI 10 ) -L 10 ) 
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G17= 1 • / ( X ( 1)>-U10) 

G18=l. /(X( ll)-lll ) 

G19= 1. / ( X( 1D-U11) 

G20= 1 . / ( X ( 12 )-L12 ) 

G2 1= 1 . /( X ( 12 ) -U1 2 ) 

G 2 2= 1 . / ( X ( 13 ) - L 1 3 ) 

G23=l . / ( X ( 13 ) -U13 ) 

S< 1» 1)=Y2*<X(4)**2 + X(9)**2)*(K11**2<-K21**2 )+2.*R*Gl**3 
S(1,2)=Y2*(X<4)*X(5)+X<9)*X(10) ) * ( K1 1 * K1 2 + K2 1 *K2 2 ) 

S( 1 * 3 ) =Y2* (X(A)*X(6)+X(9)*X(li) ) *< K1 1*K1 3+K21*K23 ) 

S( l,4)=Y2*<Kll*Hl-*-K21*H2+X(l >*X<4) *( K 1 1**2 +K2 1 **2 ) ) 

S ( 1, 5)=Y2*X<2) *X< 4 >*<K11*K12+K21*K22) 

S( It 6) = Y2* X( 3) *X ( 4 )*(K11*K13+K21*K23) 

S( 1, 7>=0.C 
S< 1,8)=0.0 

S( 1,9) =Y2*(K11*H6+K21*H7+X(1)*X(9)*(K11**2+K 21**2 ) ) 

S( 1, 1C ) = Y2 *( X( 2)*X( 9) *(K11*K12+K21*K22) I 
S( 1, 11 )=Y2*( X( 3)*X (9) *(K11*K13+K21*K2 3 ) ) 

S( 1,12 )=C. 0 
S( 1, 13 ) = C.C 

S(2,2)=Y2* ( X ( 5 )**2+X( 10) **2)*< K12**2+K 22**2 )+2.*R* 
$G2**3 

S(2»3)=Y2*(X(5)*X(6)+X( 10 )*X< 1 1 ) ) *( K 12*K 13+K22*K23 ) 

S( 2, A) =Y2*X< 1 ) *X< 5 )*(K11*K12+K21*K22 ) 

S< 2,5)=Y2*(K12*HH-K22*H2 + X< 2)*X( 5) *<K12**2+K22**2 ) ) 
S(2»6)=Y2*X(3) *X (5 )*(K12*K13+K22*K23) 

S( 2,7) =C.C 
S(2,8)=0.C 

S( 2,9) =Y2*X( 1 ) *X( 10)* ( K11*K12+K21*K22 ) 

S< 2, 10 ) = Y2 *( K 1 2*H6 + K22*H7 + X( 2 )*X( 1^) * ( K1 2**2+K 22**2 ) ) 
5(2,11 )- Y 2*X ( 3 ) *X ( 1C)*(K12*K13+K22*K23) 

S ( 2 , 12 ) =0 . 0 
S( 2, 13 )=C . C 

S( 3,3 )=Y2* ( X ( 6 )**2+X( 1 1 ) **2 ) * < K 1 3**2+K 23** 2 ) + 2 . *R* 
$G3**3 

S< 3,4)=Y2*X( 1) *X(6)*(K11*K13+K21*K23) 

S( 3,5)=Y2*X(2) *X (6 ) * ( K 1 2*K 1 3 *K22*K23 ) 

S( 3,6)=Y2*(K13*HH-K23*H2 + X(3)*X(6)*(K13**2+K23**2 ) ) 
S(3,7)=C.C 
S ( 3 , 8 ) =C . C 

S( 3, 9)=Y2*X( 1) *X( 11)*< K11*K13+K21*K23) 

S( 3, 1C ) = Y2 *X ( 2 ) *X ( 11)*(K12*K13+K22*K23) 

S( 3, 11 )=Y2*( K13*H6+K23*H7+X( 3 ) *X ( 11)* (K13**2+K 23**2) ) 

S ( 3, 12 ) = C * C 
S ( 3 * 1 3 ) =C • C 

S(4,4)=Y2*((X(1)**2)*( Kll**2+K21**2 )+K31**2 )+2 .*R* 

$< G4**3-G5**3 ) 

S(4,5)=Y2*X< 1) *X(2 )*<K11*K12+K21*K22) 

S( 4, 6)=Y2*X( 1 ) *X< 3)*( K11*K13 + K21*K23) 

S(4,7)=Y2*K31*K34 
S( 4,8)=Y2*K31*K35 
S(4,9)=C.O 
S ( A, 10 ) = C . C 
S( A, 11 )=C. C 
S( A, 12 ) = C.O 
S(A,13)=C.C 

S( 5, 5 ) =Y2* ( ( X< 2)**2) * <K12**2 + K22**2 ) +K42 **2 ) +2 . *R* 
$<G6**3-G7**3 ) 

S( 5,6)=Y2*X(2) *X(3 )*(K12*K13+K22*K23) 

S( 5,7)=Y2*KA2*KAA 
S ( 5, 8 )=Y2*KA2*KA5 
S( 5,9)=0.C 
S( 5, 10 > = 0.0 
S ( 5 , 11 )=C.C 
S(5,12)=0.0 
S(5,13)=C.C 

S( 6,6) =Y2* ( ( X (3)**2 ) * ( K1 3**2 +K23** 2 ) +K 53** 2 ) +2 . *R* 

$ ( G 8** 3-G9* *3 > 

S(6»7)=Y2*K53*K54 
S<6,8)=Y2*K53*K55 
S( 6,9)=0.0 
S ( 6 , 1C )=0.C 
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S ( 6 , 1 1 )=0.0 
S( 6, 12 )=C.C 
S ( 6 » 13 )=C.C 

S ( 7»7)=Y2*(K34**2* K44**2*K54**2 ) +2.*R*(G10**3-G11**3 ) 
S< 7,8 )=Y2*(K34*K35*K44*K45+K54*K55) 

S( 7,9) =0.0 
S( 7, 10) = C.C 
S ( 7, 1 1 )=C.O 
S(7,12)=r.C 
S(7,13) = C.C 

S< 8,8 )=Y2* (K35**2 + K45**2+K55**2) + 2.*R*(G1 2** 3-G13**3) 
S( 8, 9) =0 .0 
S( 8, 1C )=C. 0 

s ( 8 , 1 1 )=r.c 

S( 8, 12 )=C. 0 
S ( 8 , 13 )=0 . 0 

$(9,9)=Y2*((X(1 )**2)*(K11**2+K21**2)*K31**2)+2.*R* 

$ ( G 14** 3-G1 5**3 ) 



S(9,1G )=Y2*X( 1 )*X( 2)*(K11*K12+K21*K22) 

S(9,l 1 )=Y2*X{ 1 )*X (3)*( K11*K13+K21*K23) 

S( 9, 12 )=Y2*K31*K34 
S ( 9 , 13 ) = Y2*K31 *K35 

S( 10,10) =Y2*( ( X( 2)**2)*(K12**2+K22**2 ) +K42 **2 ) *2 . *R* 
$(G16**3-G17**3) 

S< 1C,1 1)=Y2*X( 2)*X(3)*(K12*K13+K22*K23) 

S( 10,12)=Y2*K42*K44 
S(1C, 1 3 )=Y 2*K 42*K4 5 

S( 11 ,1 1)=Y2*( ( X(3)**2 )*(K13**2*K23**2) +K53**2)+2.*R* 

$ ( G 18**3-G1 9**3) 

S( 11 , 1 2) = Y 2*K 53*K 54 
S( 11 ,1 3)=Y2*K53*K55 

S< 12,12)=Y2*(K 34** 2 + K44**2 + K 54**2 )+2.*R*(G2°**3-G21**3 

$ ) 

S( 12 , 1 3)=Y2* (K34*K35*K44*K45+K54*K55 ) 

S( 13, 1 3)=Y2*(K35**2+K45**2+K55**2) +2. *R* ( G22**3-G23**3 

$) 

on l i=i , nv 

DO 1 J = 1 »N V 
S(J, I ) =S ( I , J ) 

1 CONTINUE 
RETURN 
END 



SUBROUTINE MI NV( A , N,D, L . M) 

DIMENS ION A< 1 ) ,L ( 1 ) ,M( 1 ) 

DOUBLE PRECISION A , 6 , B I GA , HOLD 

PURPOSE 

INVERT A MATRIX 
USAGE 

CALL M I NV (A»N»D,L»M) 

DESCRIPTION OF PARAMETERS 

A - INPUT MATRIX, DESTROYED IN COMPUTATION AND 
REPLACED BY RESULTANT INVERSE. 

N - ORDER OF MATRIX A 
D - RESULTANT DETERMINANT 
L - WORK VECTOR OF LENGTH N 
M - WORK VECTOR OF LENGTH N 

REMARKS 

MATRIX A MUST BE A GENERAL MATRIX 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
NONE 

METHOD 

THE STANDARD GAUSS-JORDAN METHOD IS USED. THE 
DETERMINANT IS ALSO CALCULATED. A DETERMINANT OF 
ZERO INDICATES THAT THE MATRIX IS SINGULAR. 
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SEARCH FOR LARGEST ELEMENT 

0 = 1.0 

NK=-N 

00 80 K= 1 » N 
NK = NK+ N 
L ( K I =K 
M ( K ) =K 
KK=NK+ K 
B I GA = A (KK > 

DO 20 J=K , N 
IZ=N*( J-l) 

00 20 I = K » N 
I J = I Z+ I 

10 I F { DABS(B IGA)-0ABS< A( IJ )>> 15,20,20 
15 B I GA=A ( I J ) 

L ( K ) = I 
M ( K I = J 
20 CONTINUE 

INTERCHANGE ROWS 

J=L(K) 

IF(J-K) 35,35,25 
25 K I =K-N 

DO 30 1=1, N 
K I =K I + N 
H0L0 = -A( KI ) 

J I =K I-K+J 
A ( K I )= A( J I ) 

30 A ( JI ) =H0L 0 

INTERCHANGE COLUMNS 

35 I =M( K I 

IF(I-K) 45 ,45, 38 
38 JP=N* ( I— 1 ) 

00 40 J= 1 , N 
JK=NK+ J 
JI=JP+J 
H0LD=- A( JK ) 

A{ JK ) = A( JI I 
40 A( JI ) =H0LD 

DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT 
ELEMENT IS CONTAINED IN BIGA) 

45 IF(BIGA) 48,46,48 

46 D=0.C 
RETURN 

48 DO 55 1=1, N 

IF (I-K ) 5C ,55,50 
50 I K=NK+ I 

A ( IK ) = A ( IK ) / ( -BIGA I 
55 CONTINUE 

REDUCE MATRIX 

DO 65 1=1, N 
IK=NK+ I 
HOLD= A ( IK I 
I J=I-N 
DO 65 J= 1 , N 
I J=I J+N 

IF(I-K) 6C ,65,60 
60 IF(J-K) 62,65,62 
62 KJ = I J— I+K 

A ( I J ) =HOLD*A (K J ) + A ( I J ) 

65 CONTINUE 
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DIVIDE ROW BY PIVOT 



K J =K-N 

DO 75 J=1,N 

KJ=KJ+N 

IF (J-K ) 7C ,75, 7C 
70 A(KJ)=A(KJ)/BIGA 
75 CONTINUE 

PRODUCT OF PIVOTS 

0= C*BI GA 

REPLACF PIVOT BY RECIPROCAL 

A(KK )=1. C/BI GA 
80 CONTINUE 

FINAL ROW AND COLUMN INTERCHANGE 



100 

105 

108 



110 

120 



125 



130 

150 



K=N 

K= ( K- 1 ) 

I F ( K ) 150, 150, 105 

I =L ( K ) 

IF(I-K) 12C, 120,108 
JQ=N* ( K-l ) 

JR=N*< 1-1 ) 

DO 110 J = 1,N 
JK=JQ+ J 
H0LD= A ( JK ) 

J I = JR+ J 
A ( JK ) =-A ( J I ) 

A ( J I ) =H0L D 
J=M( K ) 

IF(J-K) 100,100,125 
K I =K-N 

DO 13C 1=1, N 

K I =K 1+ N 
HOLD= A (K I ) 

JI=KI-K+J 
A ( KI )=-AC J I I 
A ( J I ) =HOL D 
GO TO ICC 
RETURN 
END 
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